Astronomy & Astrophysics manuscript no. glsp 


February 2, 2008 


(DOI: will be inserted by hand later) 





Gauss-Legendre Sky Pixelization (GLESP) for CMB maps 

A.G. Doroshkevich 1 ' 2 , P.D. Naselsky 1,3,4 , O.V. Verkhodanov 1,5 , D.I. Novikov 2,6 , V.I. Turchaninov 7 , 
I.D.Novikov 1 ' 2 ' 3,8 , P.R. Christensen 1 ' 3 , L.-Y. Chiang 1,3 



in 
o 
o 



> 

m 
in 
in 
o 
m 
o 

Or 

6 

in 
c3 



Theoretical Astrophysics Center, Juliane Maries Vej 30, DK-2100, Copenhagen, Denmark 

2 Astro Space Center of Lebedev Physical Institute, Profsoyuznaya 84/32, Moscow, Russia 

3 Niels Bohr Institute, Blegdamsvej 17, DK-2100 Copenhagen, Denmark 

4 Rostov State University, Zorge 5, Rostov-Don, 344090, Russia 

5 Special Astrophysical Observatory, Nizhnij Arkhyz, Karachaj-Cherkesia, 369167, Russia 

6 Institute for Teoretisk Astrofysikk, Universitetet i Oslo, Postboks 1029 Blindern 0315, Oslo, Norway 

7 Keldysh Institute of Applied Math, Russian Academy of Science, 125047 Moscow, Russia 

8 University Observatory, Juliane Maries Vej 30, DK-2100, Copenhagen, Denmark 

Received / Accepted 

Abstract. A new scheme of sky pixelization is developed for CMB maps. The scheme is based on the Gauss- 
Legendre polynomials zeros and allows one to create strict orthogonal expansion of the map. A corresponding 
code has been implemented and comparison with other methods has been done. 
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1. Introduction 

Starting from the COBE experiment using the so called 
Quadrilateralized Sky Cube Projection (see Chan and 
O'Neill 1 1976, O'Neill and Laubscher 2 1976, Greisen and 
Calabretta 3 1993), the problem of the whole sky CMB pix- 
elization has attracted great interest. At least three meth- 
ods of the CMB celestial sphere pixelization have been 
proposed and implemented after the COBE pixelization 
scheme: the Icosahedron pixelizing by Tegmark 4 (1996), 
the Igloo pixelization by Crittenden and Turok 5 (1998, 
hereafter CT98) and the HEALPix 1 method by Gorski et 
al. 6 (1999) (with the last modification in 2003). Two im- 
portant questions mentioned already by Tegmark (1996) 
are now under discussion: a) what is the optimal method 
for the choice of the N P i X positions of pixel centers, shapes 
and sizes to provide (as good as possible) the compact uni- 
form coverage of the sky by pixels with equal areas, and 
b) what is the best way to approximate any convolutions 
of the maps by sums using pixels ? 

All the above mentioned pixelization schemes were de- 
voted to solve the first problem as accurate as possible, 
and the answer to the second question usually follows for 
the chosen pixelization scheme. 

In this paper we change the focus of the problem to 
processing on the sphere and then determine the scheme 
of pixelization. We would like to remind that pixelization 
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of the CMB data on the sphere is only some part of the 
general problem, which is the determination of the co- 
efficients of the spherical harmonic decomposition of the 
CMB signal for both anisotropy and polarization. These 
coefficients, which we call ai m , are used in subsequent 
steps in the analysis of the measured signal, and in par- 
ticular, in the determination of the power spectra, Ci, of 
the anisotropy and polarization (see review in Hivon at 
al, 7 2002), in some special methods for components sepa- 
ration (Stolyarov et al. 8 2002; Naselsky et al. 9 2003a) and 
phase statistics (Chiang et al. 10 2003; Naselsky et al. 11,12 
2003b, 2004; Coles et al. 13 2004). 

Here we propose a specific method to calculate the co- 
efficients ai m . It is based on the so called Gaussian quadra- 
tures and is presented in Sec. 2 . In this specific pixelization 
scheme correspond the position of pixel centers along the 
^-coordinate to so-called the Gauss-Legendre quadrature 
zeros and it will be shown (Sec. 5) that this method in- 
creases the accuracy of calculations essentially. 

Thus, the method of calculation of the coefficients 
ai m dictates the method of the pixelization. We call our 
method GLESP, the Gauss-Legendre Sky Pixelization. 
We have developed a special code for the GLESP approach 
and a package of codes which are necessary for the whole 
investigation of the CMB data including the determina- 
tion of anisotropy and polarization power spectra, Ci, the 
Minkowski functionals and other statistics. 
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This paper is devoted to description of the main idea 
of the GLESP method, the estimation of the accuracy of 
the different steps and of the final results, the description 
of the GLESP code and its testing. We do not discuss 
the problem of integration over a finite pixel size for the 
time ordered data in this paper. The simplest scheme of 
integration over pixel area is to use equivalent weight rel- 
atively to the center of the pixel. The GLESP code uses 
this method as HEALPix and Igloo do. 

In forthcoming papers we shall discuss our GLESP 
code extension on the processing of CMB polarization 
data, the Minkowski functionals and the peak statistics 
of CMB maps. 

2. Main ideas and basic relations 

The standard decomposition of the measured temperature 
variations on the sky, AT \9, 0), in spherical harmonics is 

oc ra—i 

AT(e,<f>)=J2 a e™Y<>m(e,(l>) (1) 

1=2 m=-l 




x 



Fig. 1. Gauss-Legendre weighting coefficients (wj) versus 
Legendre polynomial zeros [xj = cos9j) being centers of 
rings used in GLESP for the case of N = 31. Positions of 
zeros are plotted by vertical lines. 
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where Pp(x) are the associated Legendre polynomials. 
For a continuous AT(x, <f>) function, the coefficients of de- 
composition, aim, are 



1 /-27T 

dx dcj>AT(x,cf>)Y; m (x, 
1 Jo 



(3) 



where Yg denotes complex conjugation of Ye m . For nu- 
merical evaluation of the integral Eq© we will use the 
Gaussian quadratures, a method which was proposed by 
Gauss in 1814, and developed later by Christoffcl in 1877. 
As the integral over x in Eq© is an integral over a poly- 
nomial of x we may use the following equality (Press et 
al. 14 1992) 



N 



dxAT{xA)YL{x^) 
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m*m^A) (4) 



where Wj is a proper Gaussian quadrature weighting 
function. Here the weighting function Wj = w{xj) and 
AT(xj, 4>)Y(* m (xj, <p) are taken at points Xj which are the 
net of roots of the Legendre polynomial 



Pn{x 3 ) = 0, 



(5) 



where N is the maximal rank of the polynomial under con- 



sideration. It is well known that the equation Pn(xj 







has TV number of zeros in interval — 1 < x < 1. For the 
Gaussian-Legendre method EqQ, the weighting coeffi- 
cients are (Press et al. 14 1992) 



(6) 




Fig. 2. Ring center position (xj — cos8j) vs ring number 
for 2 pixelization schemes, HEALPix (solid) and GLESP 
(dashed). Figure demonstrates the case of N = 31. 



where ' denotes a derivative. They can be calculated to- 
gether with the set of Xj with the 'gauleg' code (Press et 
al. 14 1992, Sec. 4.5). 

In the GLESP approach are the trapezoidal pixels bor- 
dered by 9 and <p coordinate lines with the pixel centers 
(in the 9 direction) situated at points with Xj — cos9j. 
Thus, the interval — 1 < x < 1 is covered by N rings of 
the pixels (details are given in Sec. 3). The angular resolu- 
tion achieved in the measurement of the CMB data deter- 
mines the upper limit of summation in Eq. JQ), i < £ ma x- 
To avoid the Nyquist restrictions we use a number of pixel 
rings, ./V > 2£ max . In order to make the pixels in the equa- 
torial ring (along the 4> coordinate) nearly squared, the 
number of pixels in this direction should be N™ ax « 2N. 

The number of pixels in other rings, Ni, must be deter- 
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Fig. 3. Schematic representation of 2 types of pixelization on sphere: HEALPix (top) and GLESP (bottom). Various 
color of pixels is used to show their shape. 



mined from the condition of making the pixel sizes as equal 
as possible with the equatorial ring of pixels. 

Fig. d shows the weighting coefficients, Wj, and the 
position of pixel centers for the case N — 31. Fig. [21 com- 
pares some features of the pixelization schemes used in 
HEALPix and GLESP (see Sec. 4). Fig. |3] compares pixel 
distribution and shapes on a sphere in the mollview pro- 
jections of HEALPix and in GLESP. 

In the definition Q are the coefficients ai m complex 
quantities while AT is real. In the GLESP code started 
from the definition Q we use the following representation 
of the AT 



1=2 

+ ^ E \- a ^ Y i 
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where 
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Thus, 
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\Re[ai m ) cos(m0) - Im(ci£ m ) sin(m0)] (9) 

where P™ (cos 9) are the well known associated Legendre 
polynomials (see Gradshteyn and Ryzhik 15 2000). In the 
GLESP code, we use normalized associated Legendre 
polynomials /J": 



<2l + l {l-m)\ 
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where x = cos 8, and 9 is the polar angle. These poly- 
nomials, /™(x), can be calculated using two well known 
recurrence relations. The first of them gives /™(a;) for a 
given m and all £ > to: 



AP - 1 



21+1 (£-l) 2 - to 5 



21-3 £ 2 
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This relation starts with 
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The second recurrence relation gives /™(x) f° r a gi ven £ 
and all m < I: 



y/(i-m)(e + m+l)fF(x) = , (12) 

This relation is started with the same ff(x) and f°(x) 
which must be found with ifTTjl. 

As discussed in Press et al. 14 (1992, Sec. 5.5), the first 
recurrence relation I jlljl is formally unstable if the num- 
ber of iteration is going to infinity. Unfortunately, there 
are no theoretical recommendations what the maximum 
iteration one can use in the quasi-stability area. However, 
it can be used because we are interested in the so called 
dominant solution (Press et al. 14 1992, Sec. 5.5), which is 
approximately stable. The second recurrence relation (|12fl 
is stable for all I and m. 

3. Properties of GLESP 

Following the previous discussion we define the new pix- 
elization scheme GLESP as follows: 

— In the polar direction x = cos 9, we define Xj,j = 
1, 2, . . . , N, as the net of roots of Eq. 

— Each root Xj determines the position of a ring with 
pixel centers with ^-coordinates <pi- 

— All the pixels have nearly equal area. 

— Each pixel has weight Wj (see Eq (JBJ). 

In our numerical code which realizes the GLESP pix- 
elization scheme we use the following conditions. 

— Borders of all pixels are along the coordinate lines of 
9 and <j>. Thus with a reasonable accuracy they are 
trapezoidal. 

— The number of pixels along the azimuthal direction 
4> depends on the ring number. The code allows to 
choose an arbitrary number of these pixels. Number of 
pixels depends on the £ m ax accepted for the CMB data 
reduction. 

— To satisfy the Nyquist's theorem, the number N of 
the ring along the x = cos 9 axis must be taken as 

N ^ 2£ max ""T" L 

— To make equatorial pixels roughly square, the num- 
ber of pixels along the azimuthal axis, (f>, is taken as 
N rnax = ^(^tt / d9 k + 0.5), where k = wt(N + l)/2, 
and d6 k = 0.5(9 k+1 -e k - 1 ). 

— The nominal size of each pixel is defined as S P i xe i = 
d9k'xd(f>, where d9k is the value on the equatorial ring 
and dip = 2ir/N™ x on equator. 

— The number of pixels in the j th ring Hit/ X — X j IS 

calculated as = int(27r^/l — x 2 /S P i xe i + 0.5); 

— Polar pixels are triangular. 
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Fig. 4. Pixel size/equator pixel area vs ring- number for 
GLESP for number of rings N = 300 and N = 5000. 

— Because the number differs from 2 k where k is 
integer, we use for the Fast Fourier transform along 
the azimuthal direction the FFTW code (Frigo and 
Johnson 16 1997). This code permits one to use not 
only 2™-approach, but other base-numbers too, and 
provide even faster speed. 

With this scheme, the pixel sizes are equal inside each 
ring, and with a maximum deviation between the different 
rings of ~1.5% close to the poles (FigQJ. Increasing res- 
olution decreases an absolute error of an area due to the 
in-equivalcnce of polar and equator pixels proportionally 
to N- 2 . 

Fig.0shows that this pixelization scheme for high reso- 
lution maps (e.g. £ ma x > 500) produces nearly equal thick- 
ness d9 for the most rings. 

GLESP has not the hierarchical structure, but the 
problem of the closest pixel selection is on the software 
level. Despite GLESP is close to the Igloo pixelization 
scheme in the azimuthal approach, there is a difference 
between the two schemes in connection with the 0-angle 
(latitude) pixel step selection. Therefore, we can not unify 
these two pixelizations. The Igloo scheme applied to the 
GLESP latitude step will give too different pixel areas. 
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Fig. 5. Pixel size along polar angle (£ r , 



250) 



achieved in any pixelization scheme. To estimate this pre- 
cision we can use the expansion (CT98) 



AT" 



map 



= £ S P AT P W P (9, 4>) = £ aZ P Y em (9, 0(16) 



Lm 



> lm - , <mAT™ a *>(e,ct>)Y; mV .. 



S p AT p W;{£,m), 



(17) 



where S p is the area of the p-th. pixel. These relations gen- 
eralize Eq. J5J, taking properties of the window function 
into account. The GLESP scheme uses the properties of 
Gauss-Legendre integration in the polar direction while 
azimuthal pixelization for each ring is similar to the Igloo 
scheme and we get (see Eq(4)). 



The pixels will be neither equally spaced in latitude, nor 
of uniform area, like Igloo requires. 

4. GLESP pixel window function 

For application of the GLESP scheme, we have to take 
into account the influence of the pixel size, shape and its 
location on the sphere on the signal in the pixel and its 
contribution to the power spectrum C(£). The tempera- 
ture in a pixel is (Gorski et al. 6 1998; CT98) 



AT,, 



W p (0,(/))AT(e,(f>)dn 



(13) 



where W v 



pyw,^) is the window function for the p-th pixel 
with the area A£l p . For the window function W p (9, <f) — 1 
inside the pixel and W p (9, <fi) = outside (Gorski et al. 6 
1998), we have from Eq{IJ and Eq ((T3j> : 



AT p = J2aemW p (£,m) 



where 

W p (£,m) 

and 



dnw p {d,<t>)Y h 



(14) 



The corresponding correlation function (CT98) for the 
pixelized signal is 



(AT p AT q ) = J2C(£)W p (£,m)W;(£,m) 



(15) 



i. in 



4.1. Accuracy of the window function estimation 

The discreetness of the pixelized map determines the prop- 
erties of the signal for any pixels and restricts the precision 



W p (£,m) 



■ exp 



'2-kAx v * V 



x p -\-0.5Ax p 



in (ttto/JV*) 

(itm/N^j Jx p -0.5Ax 



dxfF(x) (18) 



where Ax p — (x p +i — x p -\)/2 with x p the p-th Gauss- 
Legendre knot and the number of pixels in the az- 
imuthal direction. This integral can be rewritten as fol- 
lows: 



x p +0.5Ax p 



dxfl n {x) 



x p -0.5Ax, 

oo 

E ^7iTT^r/ W i'(*p) 



l + (-l) fe 



k=0 



(fc + 1)! 



A.! 



fe+i 



(19) 



where f^ k \ (x p ) denotes the k-th derivatives at x = x p . 
So, for Ax p < 1 we get the expansion of 118(1: 



W^(£,m) = W^{£ 
where 

W W(£,m)~^Lexp( 



1 



/f m (A^) 2 
24/™ 



sin ^irm/N^j 



(20) 



/rwm) 



where is independent of Ax p . For the accuracy of this 
estimate we get 

SW p (e,m) = wj > 2) (£,m)~W^\£,m) (f")™(Ax p f 
W p {£,m) ~ W^\l,m) ~ 24/r 

According to the last modification of the HEALPix, 
an accuracy of the window function reproduction is about 
10 -3 . To obtain the same accuracy for the W p (£,m), we 
need to have 



(22) 



Ax p < 0.15 



£7) 

J e 



in 



(23) 
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Using the approximate link between Legendre and Bessel 5.1. Test and precision of the GLESP code 
functions for large I (Gradshteyn and Ryzhik 15 2000) 



/I™ K Jm(^ x ) 
we get: 



Ax p < Q.lbxp) ' \J m(m + 1) 

and for Ax p ~ tt/N we have from Ea.Q24|) 



8W p {i,m) >1Q _ 2 



W p (e,m) 



N 



(24) 



(25) 



For example, for N = 2t max , we obtain 
SWp(£ 1 m)/W p {£,m) ~ 2.3 x 10~ 3 , what is a quite 
reasonable accuracy for i max ~ 3000-6000. 



5. Structure of the GLESP code 

The code is developed in two levels of organization. The 
first one, which unifies F77 FORTRAN and C functions, 
subroutines and wrappers for C routines to be used for 
FORTRAN calls, consists of the main procedures 'signal' 
which transforms given values of ai m to a map, 'aim' which 
transforms a map to ai rn , 'cl2alm' which creates a sample 
of ai m coefficients for a given Ci and 'alm2cl' which calcu- 
lates Ce for ai m . Procedures for code testing, parameters 
control Kolmogorov-Smirnov analysis for Gaussianity of 
ai m and homogeneity of phase distribution, and others, 
are also included. Operation of these routines is based on 
a block of procedures calculating the Gauss-Legendre pix- 
elization for a given resolution parameter, transformation 
of angles to pixel numbers and back. 

The second level of the package contains the programs 
which are convenient for the utilization of the first level 
routines. In addition to the straight use of the already 
mentioned four main procedures, they also provide means 
to calculate map patterns generated by the Y 2 o, ^21 and 
Y22 spherical functions, to compare two sets of ai m coef- 
ficients, to convert a GLESP map to a HEALPix map, to 
convert a HEALPix map, or other maps, to a GLESP map 

Fig. [S] outlines the GLESP package. The circle defines 
the zone of the GLESP influence based on the pixelization 
library. It can include several subroutines and operating 
programs. The basic program 'cl_map' of the second level, 
shown as a big rectangle, interacts with the first level sub- 
routines. These subroutines are shown by small rectangles 
and call external libraries for the Fourier transform and 
Legendre polynomial calculations. The package reads and 
writes data both in ASCII table and FITS formats. More 
than 10 programs of the GLESP package operate in the 
GLESP zone. 

The present development of the package has also par- 
allel calculation implementation. Visualization procedures 
in OPENGL have been developed at IaO, Cambridge. 



Three tests allow us to check the code. The first of them 
is from the analytical maps 



Y 2i0 = x /—(3x 2 -l), 



*2,1 



Y-2 



2,-1 




Y? 



Yo,_ 5 



15 



-^(l-x 2 )cos(2#), 
oZir 



15 
32^ 



(1 - x 2 ) sin(2<£) 



to calculate ae m . The code reproduces the theoretical ai m 
better than 10~ 7 . 

The second test is to reproduce an analytical map 
AT(x, 4>) — Yim(x, <j>) from a given ai m These tests check 
the calculations of the map and spherical coefficients in- 
dependently. 

The third test is the reconstruction of ai m after the 
calculations of the map, AT(x,(f>) ) and back. This test 
allows one to check orthogonality. If the transformation is 
based on really orthogonal functions it has to return after 
forward and backward calculation the same ai m values. 

Precision of the code can be estimated by introduction 
of a set of aim = 1 and reconstruction of them. This test 
showed that using relation (|llfl we can reconstruct the 
introduced ai m with the precision ~ 10~ 7 limited only by 
single precision of float point data recording and with the 
precision ~ 10~ 5 for relation (|12|l . 

Fig. [7| demonstrates the accuracy of Cg calculations 
using HEALPix and GLESP 2 . 

It should be noted that unlike the HEALPix code, the 
GLESP method does not needed any iteration for calcu- 
lation of the aim coefficients and therefore is much faster. 
Our definition of the ai m coefficients is exactly the same 
as in HEALPix as an estimator of the anisotropy power 
spectrum: 



C(£) 



1 



21+1 



No I 



(26) 



5.2. Re-pixelization 



Any re-pixclization procedure will cause loss of informa- 
tion and thereby introduce uncertainties and errors. The 
GLESP code has procedures for map re-pixelization based 
on two different methods in the AT (6, 0)-domain: the first 
one consists in averaging input values in the corresponding 



2 Calculations were carried out and Fig[7]was produced by 
Vlad Stolyarov at IoA, Cambridge 
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Fig. 6. Structure of the GLESP package. 
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Fig. 7. Comparison of calculation accuracy in HEALPix 
(4 iterations) of the version 1.20 and in GLESP (no it- 
erations) methods. Number of pixels is approximately 
the same (~ 6 x 10 7 ) and calculation time is propor- 
tional to the number of iterations. Note that the com- 
parison is based on iterations of constant signal. For real 
signal (where a( rn ^ constant), the accuracy is propor- 
tional to the square root of those shown in the figure, i.e. 
ACt/d < 1(T 8 for GLESP at I < 3000 whereas - 10~ 5 
for HEALPix at all multipoles. 



pixel, the second one is connected with spline interpola- 
tion inside the pixel grid. 

In the first method, we consider input pixels which fell 
in our pixel with values AT(0i, <f>i) to be averaged with 
a weighting function. The realized weighting function is 
a function of simple averaging with equal weights. This 
method is widely used in appropriation of a given values 
to the corresponding pixel number. 




Fig. 8. Power spectra calculated for the initial HEALPix 
map (black curve) with £ m ax = 1000, N S ide — 1024, pixel 
size = 11.8026L7J', and N tot =l2 582 912, and for resulting 
re-pixelized GLESP map (red curve) with the closest pos- 
sible pixel size = 11.8038L7J', N tot =l2 581 579. Deviations 
of the power spectra at high £ illustrate the ratio of the 
HEALPix and GLESP window functions. 



In the second method of re-pixclization, we use a 
spline interpolation approach. If we have a map AT(9i, fa) 
recorded in the knots different from the Gauss-Legendre 
grid, it is possible to repixelize it to our grid AT(6' i ,fa i ) 
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using approximately the same number of pixels and the 
standard interpolation scheme based on the cubic spline 
approach for the map re-pixelization. This approach is 
sufficiently fast because the spline is calculated once for 
one vector of the tabulated data (e.g. in one ring), and 
values of interpolated function for any input argument 
are obtained by one call of separate routine (see routines 
'spline' to calculate second derivatives of interpolating 
function and 'splint' to return a cubic spline interpo- 
lated value in Press et al. 14 (1992)). 

Our spline interpolation consists of the three steps: 

— We set equidistant knots by the </>-axis to reproduce 
equidistant grid; 

— We change grid by x — cos 0-axis to the required 
GLESP grid, 

— after that, we recalculate 0-knots to the rings corre- 
sponding to the GLESP rr-points. 

Fig. |H1 demonstrates the deviation of accuracy of 
the power spectrum in a case of re-pixelization from a 
HEALPix map to a GLESP map with the same resolution. 
As one can see, for the range I < £ m ax/2, re-pixelization 
reproduces correctly all properties of the power spectra. 
For I > lmax/2 some additional investigations needs to be 
done to take into account the pixel-window function. This 
work is in progress. 

6. Resume 

We suggest a new scheme GLESP for sky pixeliza- 
tion based on the Gauss-Legendre quadrature zeros. It 
has strict expansion by the orthogonal functions which 
gives accuracy for af m -coefficients calculations below 
10~ 7 without any iterations. We realized two approaches 
for Legendre polynomials calculation using L- and M- 
methods of calculation schemes. 

Among the main advantages of this scheme are 

— a high accuracy in calculation a£ m , 

— a high speed because of no iterations, 

— an optimal selection of resolution for a given beam size, 
which means an optimal number of pixels and a pixel 
size. 

A corresponding code has been designed in 
FORTRAN 77 and C languages for procedures of 
the CMB sky map analysis. 

The a£ m calculation is the main goal. ag m -s are used 
in component separation methods and tests for non- 
Gaussianity (Chiang et al. 10 2003; Naselsky et al. 1142 
2003b, 2004). It is oriented on the fast and accurate cal- 
culation of the agm for the given resolution specified by 
the beam size. Using accurately calculated a^ m -s, one can 
reproduce any pixelization scheme by the given pixel cen- 
ters: GLESP, HEALPix, Igloo or Icosahedron. 
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